// 2009-08-15 dgl created
#include <stdlib.h>
#include "coreintegr.h"
#include <math.h>

//TOdeSys
int OdeSysInit(TOdeSys* system, int dim, TOdeFunc func,
               double* state,double time, double* params){
  if (system){
    int i;  
	  system->dim = dim;
    system->state = (double*) malloc(dim*sizeof(double));
    system->state2 = (double*) malloc(dim*sizeof(double));
	  system->params = (double*) malloc(MAX_PARAMS_NUM*sizeof(double));
    system->time = time;
    system->func = func;
    for(i=0; i<dim; i++)
      system->state[i] = state[i];
	  if (params)
	    for(i=0; i<MAX_PARAMS_NUM; i++)
        system->params[i] = params[i];
	  return C_SUCCESS;
  }
  else
    return C_FAILURE;
}

int OdeSysFree(TOdeSys* system){
  free(system->state);
  free(system->state2);
  free(system->params);
  return C_SUCCESS;
}
//dde
int DdeSysInit(TDdeSys* system, int dim, TDdeFunc func,
               double* state, double time, double* params,
               double step, double maxdelay){
  return DdeSysInitwArea(system, dim, func, state, time, params,
               step, maxdelay, NULL, 0);
  
}


int DdeSysInitwArea(TDdeSys* system, int dim, TDdeFunc func,
               double* state, double time, double* params,
               double step, double maxdelay, int *Area, int AreaCount){
  if (system)  {
    int i;  
    
    system->dim = dim;
    system->state = (double*) malloc(dim*sizeof(double));
    system->state2 = (double*) malloc(dim*sizeof(double));
    system->args = (double **) malloc(2*sizeof(double*));
  	system->params = (double *) malloc(MAX_PARAMS_NUM*sizeof(double));
    system->time = time;
    system->func = func;
    if (state)
      for(i=0; i<dim; i++)
        system->state[i] = state[i];
    if (params)
      for(i=0; i<MAX_PARAMS_NUM; i++)
        system->params[i] = params[i];

    DdeMemInit(&system->mem, step, maxdelay, system->dim);
    
    if (Area!=NULL){
      system->AreaCount = AreaCount;
      printf("Area size = %d\n",2*AreaCount*sizeof(int));
      system->Area = (int*) malloc(2*AreaCount*sizeof(int));
      for (i=0;i<AreaCount;++i)
        system->Area[i] = Area[i];
    }else{
      system->Area = (int*) malloc(2*sizeof(int));
      system->AreaCount = 1;
      system->Area[0] = 0;
      system->Area[1] = dim-1;
    }

    return C_SUCCESS;
  }
  else
    return C_FAILURE;
}

int DdeSysFree(TDdeSys* system){
  DdeMemFree(&system->mem);
  free(system->state);
  free(system->state2);
  free(system->args);
  free(system->params);
  free(system->Area);
  return C_SUCCESS;
}

//this is one-delay history
int DdeMemInit(TDdeMemory* mem, double step, double maxdelay, int dim){
 if (mem) {
   int i;
     
   TDMemNode* mn_ptr;
   mem->length = (int)(maxdelay/step);
   mem->retard = (TDMemNode *) malloc (sizeof(TDMemNode));
   mem->retard->v = (double*) malloc (dim*sizeof(double));
   mn_ptr = mem->retard;
   for (i=1; i<mem->length; i++){
     mn_ptr->next = (TDMemNode *) malloc (sizeof(TDMemNode));
     mn_ptr = mn_ptr->next;
     mn_ptr->v = (double*) malloc (dim*sizeof(double));
   }
   mn_ptr->next = mem->retard;
     
   return C_SUCCESS;
 }
 else
   return C_FAILURE;
 
}

int DdeHistoryFillwConstant(TDdeMemory* mem, double* xinit, int dim, double step){
  if (mem)
  {
    TDMemNode* mn_ptr;
    int i,j;
    
    for (i = 0; i<dim; i++)
      mem->retard->v[i] = xinit[i];
    mem->retard->t = - step * (double)mem->length;
        
    mn_ptr = mem->retard->next;
    j = 1;
    while (mn_ptr != mem->retard){
      for (i = 0; i<dim; i++)
        mn_ptr->v[i] = xinit[i];            
      mn_ptr->t = -step * (double)(mem->length-j);
      mn_ptr = mn_ptr->next;
      j++;
    }
         
    return C_SUCCESS;
  }
  else
    return C_FAILURE;
}
    
   
int DdeMemFree(TDdeMemory* mem){
  if ((mem)&&(mem->retard))
  {
    TDMemNode* mn_ptr=mem->retard->next;
    TDMemNode* killer;
    while (mem->retard != mn_ptr){
      killer = mn_ptr;
      mn_ptr = mn_ptr->next;
      free(killer->v);
      free(killer);
    }  
    free(mn_ptr->v);
    free(mn_ptr);
  }           
  return C_SUCCESS;
}


//MethodData
int MethodDataInit(TMethodData* data, int dim, intmethod methname, double step, double atol, double rtol){
  if (data){
    int i, num;  
    
    switch (methname){
      case IM_EULER:
        num = 1;    
        break;
      case IM_RK4:
        num = 5;    
        break;
		  case IM_DOPRI54:
        num = 8;    
        break;
      case IM_RADAUIIA:
        num = 6;    
        break;

      default:
        return C_FAILURE;
    }
    data->dim = dim;
    data->num = num;
    data->step = step;
    data->stor = (double**) malloc(num*sizeof(double*));
    for (i=0; i<num; i++)
      data->stor[i] = (double*) malloc(dim*sizeof(double));
    data->state_ref = 0;
	  data->atol = atol;
	  data->rtol = rtol;
	  data->rejected = 0;
	  data->accepted = 0;
#ifdef __DGL_DEBUG
	  data->log = fopen("metod.log", "w");
#endif
    return C_SUCCESS;
  }
  else
    return C_FAILURE;
}

int MethodDataFree(TMethodData* data){ 
  int i;
  for (i=0; i<data->num; i++)
    free(data->stor[i]);
  free(data->stor);
#ifdef __DGL_DEBUG
  fclose(data->log);
#endif
  return C_SUCCESS;
}


int OdeIntegrateNstep(TOdeMethod integrator,TOdeSys* system, int step_num){
  int i;
  for(i=0;i<step_num;i++)
    integrator(system); 
  return C_SUCCESS;
}

